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Abstract 

This article deals with dynamical systems depending on a slowly varying 
parameter. We present several physical examples illustrating memory effects, 
such as metastability and hysteresis, which frequently appear in these systems. 
A mathematical theory is outlined, which allows to show existence of hysteresis 
cycles, and determine related scaling laws. 
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1 Introduction 



There exist many instances where the dynamics of a system depends on a parameter 
which varies slowly in time. This parameter is often controllable by the experimen- 
talist, who can modify it at will. A well-known example of this situation is that of 
a ferromagnet on which is imposed a low frequency magnetic field. One can also 
think of chemical reactions occurring in a reactor in which the flux of the injected 
chemical substances is varying slowly, or the Couette-Taylor experiment in hydro- 
dynamics where the speed of rotation of the inner cylinder is slowly modulated. In 
other circumstances, the parameter is not controllable, but certainly influences the 
dynamics of the system of interest. As examples of this situation we could mention 
the case of the impact of solar light on the thermal convection in the atmosphere, 
or the seasonal (or even climatic) effect on the dynamics of populations. 

One of the most interesting phenomena observed in systems with an adiabatically 
varying parameter is the familiar one of hysteresis. Recently, there has been a 
renewal of interest in this old problem, both from a theoretical and an experimental 
point of view. Several authors ||RKP| , EE] , |J GRM| , |Hhfc|| have particularly analysed 
properties of the hysteresis cycle, such as its area, which appears to scale in a 
nontrivial way with the adiabatic parameter. 

In this article, we concentrate on dynamical systems with a finite number of de- 
grees of freedom, depending on a parameter in such a way that the system undergoes 
bifurcations when the parameter is considered to be static. The static (or "frozen") 
situation corresponds to measurements made, in principle during very long times, 
at successive fixed values of the parameter. We then ask what is happening when 
the parameter is varying slowly in time, instead of being kept fixed. This question is 
closely related to the opposite one: can the static bifurcation diagram be determined 
experimentally by varying the parameter slowly in time (a possible temptation for 
the impatient experimentalist)? 

We have recently developed a coherent mathematical framework to deal with 
adiabatic systems, in particular to show existence of hysteresis cycles and determine 
their scaling laws ||Berj| . The purpose of this article is to explain these methods by 



illustrating them on a few concrete physical examples. 

We begin, in Section [| by presenting the most important features of one- 
dimensional (ID) systems, which are illustrated by a few generic examples. We 
discuss in particular a simple geometric method to determine scaling laws near bi- 
furcation points. 

In Section we use the Lorenz model to illustrate the phenomenon of bifurcation 
delay. When translated into the language of Rayleigh-Benard (RB) convection, this 
phenomenon means that the slow and periodic variation of the temperature gradient 
in time leads to a delayed appearance of convection rolls and to hysteresis. The 
Lorenz model being a good approximation close to the instability threshold, since 
it contains two dominant modes of the bifurcation, this delay should be observable 
in the real RB convection. To explain the delay, we introduce two new methods 
especially designed for n-dimensional (nD) systems: dynamic diagonalization and 
adiabatic manifolds. 
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In Section |], we present a simple mean field model for the dynamics of a ferro- 
magnet in a slowly oscillating magnetic field. In the ID case, we discuss the concept 
of dynamic phase transition introduced in [ [T0|| , and derive a scaling law for the hys- 
teresis area. In the 2D case, we examine the effect of anisotropy on the mechanism 
of magnetization reversal and the shape of hysteresis cycles. 

In Section BL we discuss a simple mechanical system (which was introduced in 
PK|| ), displaying chaotic instead of periodic hysteresis. This phenomenon depends 
only on a few qualitative features of the system, and should be observable in a larger 
class of nonlinear oscillators including inertia and involving a symmetry breaking 
bifurcation. 

Finally, Section |6] is dedicated to some examples of the effect of eigenvalue cross- 
ings. These crossing give rise to an effective interaction between otherwise indepen- 
dent modes, which is essential in the sense that it cannot be eliminated by a change 
of variables. The interaction may, however, be delayed in certain cases. 

Throughout this text, we use the following mathematical setting. The "frozen" 
dynamical system is supposed to be described by a family of ordinary differential 
equations 

H x 

— — F(x, A), xeR n , XeR. (1) 
The associated adiabatic system is given by 

^ = F(x,X(et)), (2) 

where A(r) is a given function, and e is the small adiabatic parameter. It is useful 
to introduce the slow time r = et, so that (0) can be rewritten as 

e^ = F(x,X{r)), (3) 

or, in short form, ex = f(x,r). We denote by (-|-) the usual scalar product in R n 
and by ||-|| the Euclidean norm. 

There is a large literature on singular perturbed problems of this type. Results 
on linear systems can be found in [ |Wa| | . For a review of results on dynamic bi- 



furcations, see |Pen|| . In particular, the phenomenon of bifurcation delay has been 



rigorously analysed in two important articles by Neishtadt [|Nel| , [Ne2| . Certain hys- 



teresis phenomena in slow-fast systems similar to the Van der Pol equation have 
been analysed in |[MKfe|| . Our geometric method to determine scaling exponents, 



as well as the procedures using dynamic diagonalization and adiabatic manifolds, 
however, are new to our knowledge. Here we only outline some essential ideas of 
these methods, detailed proofs can be found in |[Berj . 
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2 One— dimensional systems 



In this section, we will consider one-dimensional (ID) adiabatic equations of the 
form 

ex = F(x, A(r)) = f(x,r), x, r e R , (4) 

where the dot denotes derivation with respect to r, and f(x, r) is assumed to be an 
analytic real-valued function (weaker results as those stated below hold for differ- 
entiable functions). 

The static bifurcation diagram of (f|) is obtained by determining the solutions of 
f(x,r) = 0, which are generically curves subdividing the plane into regions where / 
is positive or negative. Let x*(r) be such an equilibrium curve. The implicit function 
theorem tells us that if the linearization a(r) = d x f(x*{r),T) does not vanish, then 
x*(t) is a smooth curve. It corresponds to stable solutions if a(r) is negative, and 
to unstable ones of a(r) is positive. 

In such a situation, one can prove the existence of a particular solution x{r) of 
the adiabatic equation @ tracking the curve x*(r) at a distance of order e: 

x{t) = x*(r) + 0(e). (5) 

Moreover, this solution admits an asymptotic power series in e, which does not 
converge in general, but admits, however, an optimal truncation at exponentially 
small order: 

N(e) 

x(t) = x*(t) + ZiOy + 0(e~ 1/C£ ), N{e) = 0(l/e). (6) 
j=i 

We call x(t) an adiabatic solution associated with the equilibrium branch x*(r). 
Other solutions of (||) are attracted or repelled exponentially fast by adiabatic ones, 
and tend to switch between the neighborhoods of different equilibrium branches in 
a time of order e|liie|. As long as there are no bifurcation points, the solutions thus 
remain most of the time close to equilibria, and there is no hysteresis in the system. 

Let us now consider the effect of bifurcations. At these special points, several 
equilibrium branches may meet, causing the solutions to choose between several pos- 
sible directions, which is the basic mechanism of hysteresis. Moreover, equilibrium 
branches are in general no longer tracked at a distance of order e, but at a distance 
scaling in some other, nontrivial way with e, which we now show how to compute. 

If the origin is a bifurcation point of /, we can write in some neighborhood 

f(x,r)= c nm x n T m , c 00 = cio = 0. (7) 

n,m ^ 

Assume that f(x, r) admits an equilibrium branch scaling as x*(t) ~ \r\ q (we use this 
notation to indicate that C-Irl 9 ^ x*(t) ^ c + |r| g , where c± are positive constants 
independent of r and e) . A standard result of bifurcation theory states that — q is 
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Figure 1. Newton's polygons for the most generic bifurcations discussed in the 
text. Dots mark points for which c nm ^ 0. The slope of segments correspond to 
the possible exponents q of equilibrium branches through the bifurcation point. 
The ordinate at 1 of these segments is the exponent p of the linearization. 



necessarily equal to the slope of a segment of Newton's polygon. This polygon 
is constructed as the convex envelope of the set of points (n, m) such that c nm ^ 
0, completed by a horizontal and a vertical (Fig. |l|). The linearization a(r) = 
d x f(x*(r),r) scales generically as \t\ p , where p is the ordinate at 1 of the tangent 
to Newton's polygon with slope —q. 

These two easily determined numbers p and q are usually sufficient to characterize 
the scaling behaviour at leading order in e. In fact, different behaviours take place 
in an inner region |r| ^ e 1 ^^ 1 and in an outer region \r\ ^ e v ' p+l . In particular, 
if x*(t) ~ \r\ q is a decreasing stable branch arriving at the bifurcation point, above 
which / is negative, one can show that 



X(T) 



x*(r) 



£ \ T \q-P-i for r ^ - e 1/p+1 , 



r?/P+l 



for < r < 0. 



Combining a local analysis around bifurcation points with a global analysis, 
which is usually easy in ID, one can determine the qualitative properties of dynamics. 
In particular, if A(r) is a periodic function, one can construct the Poincare map 
(which is necessarily a monotonous function) in order to prove existence of hysteresis 
cycles and determine their scaling laws. 

Let us illustrate this procedure on the simple model equation given by 



F(x, A, ji) = —fix — x 3 + A. 



(9) 



Mathematically, this function is a generic two-parameter perturbation of the vector 
field —a; 3 [ |HK| . Physically, it describes the overdamped motion of a particle in a 
Ginzburg-Landau potential 



$(x, A, p,) = t^x 2 + ~ 



(10) 



where fi = T — T c represents the difference between the temperature and its critical 
value, and A is an external field. The quartic potential described by the first two 
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Figure 2. (a) Orbits, in the (A, a^-plane, of the system ex = — fix — x 3 + A(r), for 
/j, = 1 and e = 10~ 3 / 2 , 10~ 2 and 10~ 5 / 2 . The solutions are attracted by a periodic 
orbit, enclosing an area of order e. (b) Same as (a), but for \x = 0. At A = 0, the 
periodic orbit is at a distance of order e 1 / 5 of the equilibrium point, and it encloses 
an area of order e 4 / 5 . 



terms is fairly generic in physical systems presenting the symmetry x i— > —x, while 
the linear term is the simplest possible asymmetric perturbation. 

We begin with the situation where \i — 1 is fixed, and A(r) = sinr is slowly 
oscillating. The equation F(x, A, 1) = admits a single, stable equilibrium branch 
x*(X), given implicitly by x*(A) 3 + x*(X) = A. All solutions are attracted by a 
periodic solution x(r) = x*(A(r)) + 0(e), enclosing an area of order e (Fig. |2]a). In 
the adiabatic limit e — > 0, this area vanishes and there is no hysteresis. 

If n = and A(r) = sinr, the unique equilibrium branch x*(X) = A 1 / 3 admits 
the origin as a bifurcation point. Using Newton's polygon or a direct calculation, we 
find that the exponents determining the scaling behaviour are q = | and p — | (Fig. 
|l]a). The orbits are attracted by a periodic one, crossing the x-axis at a distance of 
order e 1//5 from the origin, and enclosing an area 

A(e) w £ 4/5 (11) 

(Fig. |b). The cycle still collapses with x*(A(r)) in the adiabatic limit, but with 
a much slower rate. These exponents have been found in ||HL&] , |GBS|| using other 
methods (which are less general than ours). 

If fi — — 1 and A(r) = sinr, there are two bifurcation points at (±A C , =Fx c ), where 
A c = a/4/27 and x c = a/1/3. Two stable branches x^_(r) and one unstable branch 
Xq(t) meet at these bifurcation points (Fig. §a). Since they are crossed with nonzero 
velocity, Newton's polygon shows that the associated exponents are q = p = \ (Fig. 
[l|b). In fact, close to these points, the dynamics in translated coordinates is governed 
by the equation 

ey = —r — y 2 + higher order terms. (12) 

Solutions cross the y-axis at y ~ e 1 ^ 3 . Using the scaling y = £ 1//3 z and r = e 2 ^ 3 a, 
one shows that y remains of order e 1 ^ 3 until a time r* w e 2 / 3 , and then quickly leaves 
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Figure 3. (a) Same as Fig. for /j, = — 1. The periodic orbit lies at a distance 
at most 0(e 1 ^ 3 ) from a limiting hysteresis cycle, composed of stable equilibrium 
branches and two vertical lines. The enclosed area is A(e) = A(0) + 0{e 2 ^ 3 ). (b) 
Hysteresis cycle of the equation ex = [i(t)x — x 3 . 



the bifurcation region to jump on the other stable branch. The orbits are attracted 
by a hysteresis cycle with area satisfying 

A{e) - -4(0) « e 2/3 , (13) 

where A(0) — § is the area situated between £+(A) and x* (A) for — A c < A < 
A c . This time, the hysteretic behaviour persists in the adiabatic limit. The main 
contribution of order e 2 ^ 3 to the excess area comes from the delayed jump. The 
scaling law ( |T3"D was also obtained in [ JURM using an exact solution of fll^p (without 



the higher order terms). 

As a final example, let us consider the situation where A = and fi(r) is the 
varying parameter. The static bifurcation diagram displays a pitchfork bifurcation 
at the origin, involving the branches x = and x = ±a/A, for which q = | and 
p = 1 (Fig. 0c). An important new phenomenon is bifurcation delay: near x = 0, 
the dynamics is essentially governed by the linearized equation 



ex 



i r 

)jl{t)x x(r) = exp[- / fj,(s) ds\x(r ). (14) 



Starting at a time where \i < 0, x(r) remains exponentially small as long as the 
integral in (14] ) is negative, which is true for a while beyond the instant where 



fx becomes positive. If the solution finally jumps on the stable branch and \x is 
decreased again, this branch is followed adiabatically, which leads to hysteresis (Fig. 
0b). The area of the cycle follows the scaling law 

A(e) - A(0) w e 3/ \ (15) 

where A(Q) depends on the bifurcation delay time. 

These examples indicate that one-dimensional systems are relatively well under- 
stood. In fact, our method to determine scaling laws is general, and quite straightfor- 
ward to apply to other bifurcations than the special cases described here. Moreover, 
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one can expect that when a mode of a larger system undergoes bifurcation, the 
dynamics of the associated adiabatic system will be governed by an effective ID 
equation, describing the motion of this particular mode, which explains why the 
same scaling laws are observed for more complicated systems. The Lorenz model 
discussed in the next section illustrates this reduction of variables. 

3 Bifurcation delay in the Lorenz model 

Let us now turn to the behaviour of higher-dimensional systems 

ex = F(x,X(t)) = f{x,r), reR. (16) 

Adiabatic solutions with an expansion of the form (§) still exist in the vicinity of 
hyperbolic equilibria, i.e., branches of equilibrium points around which the lin- 
earization of / has no purely imaginary eigenvalues. 

The behaviour of neighbouring solutions is far more complicated to analyse than 
in the ID case. One can, however, obtain valuable informations by analysing the 
linearized system first, before dealing with nonlinear terms using appropriate tools 
such as invariant manifolds. 

We illustrate these techniques on the Lorenz model with slowly varying temper- 
ature gradient r(r): 

ex\ = a{x 2 — X\) 

6X2 = r{r)x\ — X2 — X1X3 (17) 

EX 3 = —bx 3 + XiX 2 , 

where we assume that b, o > 0. This model has been introduced as an approximation 
to Rayleigh-Benard convection, but also describes other systems such as lasers. It is 
well known that if r ^ 1 is fixed, the origin is a globally asymptotically stable fixed 
point, whereas when r > 1, the origin is hyperbolic and two new equilibria C± = 
{±y/b{r- 1), ±y/b(r-l) , r — 1) appear, which correspond physically to convection 
rolls. 

We will study this system when r(r) is slowly oscillating around r = 1, and stays 
well below the chaotic region. It can be written in compact form 

ex = A(t)x + b(x), (18) 

where b(x) is quadratic. Note that the identically zero function is a particular 
solution of this equation. The matrix A(r) has three eigenvalues 

ai, 2 (r) = -i(«7 + l)±s(r), a 3 = -b, s{r) := |a/ (<t + l) 2 + 4er(r(r) — 1). (19) 

The eigenvalues a 2 and a 3 are always negative, while ai(r) has the same sign as 
r(r) — 1. We thus expect that the motion will essentially follow the eigenspace of 
ai(r). For this reason, we will construct a change of variables isolating this particular 
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direction. To do this, we begin by searching a linear transformation which should 
diagonalize the linear part of flTSp. To this end, observe that if S(r;e) is a matrix 
satisfying 

eS = AS - SD (20) 
where D(t; e) is diagonal, then the change of variables x = Sy transforms (|T8|) into 

ey = D[r)y + S- 1 b{S V ). (21) 



The key point is that we can prove the existence of a bounded solution of (p0|). 
admitting asymptotic series 



S(t; e) = S (t) + eS'i(r) + e 2 S 2 {r) + ■ ■ 
D(t; e) = D (r) + e£>i(r) + e 2 D 2 (r) + • 



2o,_x, (22) 



which can be truncated to exponentially small order, just as the adiabatic solution 
(§]).[] In particular, Sq(t) is the matrix diagonalizing A statically, and the entries of 
Dq(t) are eigenvalues of A(r). The proof uses the fact that ai(r) ^ a 2 (r). 

In the specific case of the Lorenz equations, a linear transformation given to 
leading order by 

x 1 = a{y + z 1 ) + 0{e), x 2 = + z x ) + s(y - z x ) + 0(e), x 3 = z 2 (23) 

yields the equation 

sy = di(r)j/ + &i(!/,*,T) 
ei = £> 2 (t)^ + b 2 (y,z,r), 

where c?i(r) = cti(r) + C(e), D 2 (t) is diagonal with entries a 2 (r) + 0(e) and 03, and 
61, b 2 are quadratic. 

In order to deal with the nonlinear terms, we introduce invariant manifolds. 
Consider the partial differential equation 

ed T v(y, r) = D 2 (r)v(y, r) + b 2 (y, v, r) - d y v(y, r) [di(r)j/ + 61 (y, u, r)] . (25) 

It can be shown that this equation admits, in a neighborhood of y = 0, a solution 
v (?/> r ) — 0(?/ 2 )- When r > 1, it lies at a distance of order e from the instantaneous 
unstable manifold of the origin, but it can be continued to times where r < 1. The 
change of variables z = ( + v(y,r) transforms the second equation of ([24]) into 

eC= [D 2 (r) + f3 2 (y,C,r)]C, (26) 

where (3 2 is of order \y\ + ||£||. This equation admits ( = as invariant manifold. A 
similar change of variables y = 77 + iz(£, r) transforms the first equation into 

87]= [di(r)+#ifaC,T)]77, (27) 



1 These matrices are not unique, since every column of S 1 can be multiplied by a function of 
time, which will of course affect terms of order e in D. 
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Figure 4. Rotation frequency of convection rolls x\ as a function of the periodi- 
cally varying temperature difference r in the Lorenz model. After a first transient 
cycle, the motion settles on a hysteresis cycle, on which x\ increases rapidly and 
decreases slowly. 



defining a stable manifold separating the basins of attraction of C + and C_. 

Since D 2 has negative eigenvalues, one easily shows that £(r) goes to zero expo- 
nentially fast. Thus the effective dynamics will take place on the invariant manifold 
z = v(y, t), where it is governed by the scalar equation 

ey = d 1 (T)y + h(r)y 3 + 0(y 5 ), 

a 2 /a -I \ (28) 

(y = 1] on the unstable manifold ( = 0). As r(r) is varied periodically around r = 1, 
di(r) = ai(r) + 0(e) changes sign. The situation is thus very similar to the last 
example of Section ||. Assume that di{r) becomes negative at r , positive at n, and 
negative again at tq + 1. Asymptotically, the solution will stay close to the origin 
for r + n^r^r + n, where f > t x is the delay time denned by the relation 

rf 1 (r)dr = (29) 

(this time exists if the average of d\{j) over one period is positive, otherwise solutions 
stay indefinitely close to the origin). For f + n<T<r + n + l, the asymptotic 
solution follows C + or C_ adiabatically. Thus the bifurcation delay leads once again 
to hysteresis (Fig. ||). 

When r(r) is varied back and forth (at least when r does not become too large), 
the solution always follows the same equilibrium (which one it chooses depends 
on the initial condition). In the case of Rayleigh-Benard convection, x± measures 
the rotation frequency of convection rolls. When r(r) is increased, these rolls will 
appear suddenly, with a positive frequency, at some r c = r(f) > 1. When r(r) is 
decreased again, they slowly decelerate to disappear smoothly as r becomes smaller 
than 1. The rolls will always turn in the same direction. We believe that it would 
be interesting to try to observe this delay experimentally. 



Ti 
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4 Hysteresis in mean field ferromagnets 



Hysteresis in ferromagnets has been known and studied experimentally for a long 
time. Interest in a microscopic understanding of hysteresis and associated scaling 
laws has been renewed by the numerical study of JKKP]. The internal dynamics of 



ferromagnets, however, is so complicated that its modeling by ordinary differential 
equations is not obvious. We will consider here a simple Curie- Weiss model, which 
can be described by an effective mean field equation. 
Consider the Hamiltonian 

where A is a subset of Z d with N sites, the spins <Ji are unit vectors in R n , J is 
a fixed coupling matrix and h the magnetic field. We introduce a stochastic spin 
flip dynamics of Glauber type. The detailed balance condition |[Ka|| is satisfied by a 
transition probability by unit time of the form 

= E II 5 K - J {aiMtr)) 9(hM)l 

1 (31) 
hi(a):=h+— Jj2 a j> 

where (3 = T~ l is the inverse temperature, and g(h) an arbitrary function. 

To derive a deterministic equation of motion, we consider a sequence of systems 
with N sites, iV — > oo. Under appropriate assumptions on the initial probability 
distribution, one can derive in the thermodynamic limit a deterministic equation for 
the magnetization m of the form| 

(1777 

— = -m + P(Jm + h)F n (J3\\Jm + h\\), (32) 
at 

where F n (x) depends on the dimension n of the spins. In particular, 

th r 11 

F 1 (x) = , F 2 (x) = ---x 2 + 0(x 4 ). (33) 

x 2 lb 

It can be shown that corrections to m(t) resulting from a finite iV are of order iV -1 / 2 , 
and obey a Langevin equation | |Mal| , |Ma2 1 . 



We will now describe some interesting properties of the ID and 2D models in a 
slowly oscillating magnetic field, paying attention in particular to the mechanism of 
magnetization reversal and its influence on the shape of hysteresis cycles. 
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Figure 5. Solutions of (34) illustrating the phenomenon of "dynamic phase tran- 
sition". When the amplitude h of the magnetic field is smaller than the critical 
field h c , the magnetization oscillates around a nonzero average, and encloses an 
area of order eh (a). When h is larger than h c , the average magnetization is zero, 
and the periodic solution encloses an area A{e) = .4(0) + 0(e 2 / 3 ) (b). 

4.1 One— dimensional case: dynamic phase transition 

For ID spins, the adiabatic equation of motion can be written as0 

em = -m + th(3(Jm + h(T)), (34) 

where we will consider a periodic magnetic field of the form 

h{j) =/i sin(27rr). (35) 

For positive inverse temperature (3 = T~ x , we may rescale the variables in such a 
way that J = 1. If (3 < 1, there is no static hysteresis. If f3 > 1, the static bifurcation 
diagram is similar to the equation ex = —x + x 3 + A(r) discussed in Section it 
contains two saddle-node bifurcations (±/i c , T m c), where 

m c (T) 



h c (T) 



— T Argthm c . 



(36) 



As we have seen in Section pj, when h Q ^> h c orbits are attracted by a hysteresis cycle 
with zero average magnetization (Fig. |5|b) and area *4(0) + 0(e 2 ^ 3 ). This scaling law 
was already obtained in | JGRM |. When h Q < h c , the magnetization never sees any 
bifurcation point, and it follows asymptotically a cycle of width 0(e) with nonzero 
average magnetization (Fig. |5]a). 

This phenomenon was observed numerically and called "dynamic phase transi- 
tion" by |[rO|| , who also called "ferromagnetic" or F— region the domain of (T, ho)- 
plane where the asymptotic cycle has a nonzero magnetization, and "paramagnetic" 
or P— region the domain where it has zero average magnetization. In the adia- 
batic limit, these regions are delimited by the line h = h c {T). For positive e, there 

2 This is the simplest equation, obtained for a particular choice of g. Other choices yield a 
multiplicative factor in the right-hand-side. 

3 We neglect here terms of order e stemming from the slow time dependence of h in the derivation 
of the equation of motion. 
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Figure 6. If the amplitude of the magnetic field is equal to h = h c + 5, the 
motion near the turning point is governed by the Ricatti equation (^). If S = 0, 
it describes a transcritical bifurcation, and adiabatic solutions follow the upper 
branch (a). This means that we are still in the F-region. For positive 8, we show 
that this behaviour subsists as long as 8 = 0(e). In (b), trajectories are shown for 
two different values of e. If e < 8, the solution escapes from below after a delay of 
order e 2 / 3 ^ 1 / 6 , and we have reached the P-region. 



may be a small overlap between these regions, where a symmetric P-cycle and an 
asymmetric F-cycle coexist. 

We claim that for small e and T = f3~ l < 1, the F-region grows by a distance 
of order e, the P-region shrinks by an amount of the same order (but may overlap 
the F-region), and the corresponding cycles obey the scaling laws 

F-cycle: ho) ~ h e if h < h c , eln|e| if h = h c , 

P-cycle: A(e, h ) « A + e 2/3 (h - h c ) 1/3 . ' °' ' 

Let us indicate how we obtain these scaling laws. Assume that h = h c + 5. After 
translating the coordinates to the bifurcation point (—h c , m c ) and scaling them in 
a proper way, equation ( |3l] ) becomes 

ey = —y 2 — 5 + r 2 + higher order terms. (38) 

For 5 = 0, this equation displays a transcritical bifurcation at the origin. For positive 
<5, it splits up into two saddle-node bifurcations, with a gap of width 2^5 (Fig. g). 
If 5 < 0(e), the transformation r = \f5~o, y = y/ez yields the equation 

r] Z , 

^1 = -ez 2 + £ 3 (a 2 - 1), e= v^/i, (39) 
da 

which can be used to show that z cannot move enough to slip through the gap, so 
that we are in the P-region. 

If 5 > 0(e), the transformation r = y/5(a — 1), y = V~5z gives 

( £ r 1 )^ = - 2 2 -a + a 2 , (40) 
do" 
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Figure 7. (a) Evolution of y(r), solution of equation (f43|). The magnetization 
quickly rotates at times r nj determined recursively by the relation r n+1 = § + 
f (t„ — S). (b) The plot of mi as a function of hi shows the asymptotic hysteresis 
cycle, which is determined solely by the delay times. Thick solid lines, dashed and 
dotted lines represent respectively sinks, saddles and sources of the static system. 
Due to bifurcation delay, the magnetization follows the hyperbolic branch for some 
time, but ultimately rotates around the unstable origin. 



which is exactly the equation studied in Section FJ. In particular, the trajectory 
slips through the gap after a time delay of order ^(ed^ 1 ) 2 ^ 3 = s 2 /' 3 5~ 1 / 6 . During 
this time, the magnetic field has reached a value ho + 0(e 2 ^ 3 S 1 ^ 3 ), which implies the 
scaling relation fl37D - 

Finally, when 8 m e, the trajectory may behave in either way. A more careful 
analysis of the Poincare map shows that even though there is a small region where 
stable F- and P-cycles can coexist, the transition is sharp, in the sense that the 
average magnetization jumps discontinuously from one cycle to another. In ||TO|| , 



a smooth transition, where the magnetization goes to continuously, has been ob- 
served for larger values of e. 

4.2 Two— dimensional case: effect of anisotropy 

If we retain only the leading terms in (p2|) , we obtain the Ginzburg-Landau equation 

em = ((3J-l)m-^(3Jm\\(3Jm\\ 2 + (3h } (41) 

which describes the linearly driven, overdamped motion of a particle in a sombrero- 
shaped potential. We shall assume that h is parallel to an eigenvector of the sym- 
metric matrix J. In the subspace of this eigenvector, the equation reduces to the 
previously studied case, and a minimal field amplitude h c is necessary to reverse 
magnetization. When the magnetization has a transverse component, however, it 
can also turn around the potential maximum, for a much smaller field amplitude. 
We are thus going to focus on this situation. 

In the isotropic case, we may choose J = 11 to obtain the equation 



em 



{(3 - l)m - l(3 3 m\\m\\ 2 + ph(r), (42) 
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Figure 8. (a) Evolution of ip in the anisotropic case with 7 > 1. Due to bifurca- 
tion delay, ip spends some time near or it, even when these points are unstable. 
It always drops back, however, to the transverse stable position, (b) The resulting 
hysteresis loop looks triangular. Curved lines are longitudinal equilibria, and the 
straight line represents transverse branches. 



where we take a magnetic field h{r) = (7ii(t),0), with h\(r) = h sin(27rr). It is 
useful to write this equation in polar coordinates, with m = (r cos ip,r simp), to get 



er 



eip 



{(3 - l)r - |/3 3 r 3 + ph^r) cosf 
P 



—hi(r) sin if . 

r 



(43) 



If hi(ro) < 0, the magnetization settles near the left equilibrium, determined by 
if — 7r and r = r„(r ), the largest solution of (f3 — l)r — |/5 3 r 3 — /5/ii(r ) = 0. 
When the field becomes positive, the phenomenon of bifurcation delay causes if to 
remain for some time in unstable equilibrium near n, until it switches to at a time 
7~i = \I/(t ), defined by 



W(to) 



To 



Kr) 
r_(r) 



dr = 0. 



(44) 



Because of the symmetry, the next time of delayed magnetization reversal is then 
given by r 2 = | + ^(ti — |) (Fig. 0). Subsequent reversal times are determined by 
the recursive formula r n = | + ^f(r n — |). It turns out that this self— determined 
bifurcation delay finally settles at a fixed point of the map r 1— > ^(t) — |. 

We now turn to the anisotropic case where J — (07)- m t ne coordinates 
m = (r cos(/9,7 _1 r siny)), the second equation of (fl"3|) becomes 



e0 = /5(1 — |/3 2 r 2 )(7 — 1) sin if cosy? — — /ii(r) sin<£> 



(45) 



The case 7 < 1 is not very interesting, since the anisotropy enhances the effect of 
the magnetic field, and tends to align the magnetization with it. If 7 > 1, a new 
stable transversal equilibrium exists for small magnetic field. Its coordinates are 
determined by the relations ||/3Jm|| 2 = 2(1 — /3 _1 7 _1 ) and (1 — 7~ 1 )m 1 = fih\. The 
resulting hysteresis cycle is composed of two triangular loops (Fig. g), since after 
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Figure 9. Phase portraits of the rotating pendulum for different values of the 
rotation frequency 0. 

leaving the unstable position ip = or tt, the magnetization drops to the transverse 
branch, which it follows until merging with the longitudinal branch. 

We point out that if the magnetic field is slightly tilted with respect to the 
eigenvectors of J, the pitchfork bifurcations in Fig. |8| transform into saddle-nodes, 
which suppresses the bifurcation delay. The result is that instead of oscillating back 
and forth, the magnetization performs full circles, always in the same direction. 

5 Chaotic hysteresis of a rotating pendulum 

The examples we have considered up to now all described an overdamped, effectively 
one-dimensional motion, which displayed hysteretic, but not chaotic properties. We 
present here an example taking into account inertia, which turns out to have far 
more complicated dynamics. 

Consider a mathematical pendulum mounted on a rotating table, turning with 
angular frequency Q. The pendulum is subject to weight, friction and a centrifugal 
torque, so that its equation of motion can be written in dimensionless variables 

Q = V , s 

„ (46) 
p = —2"fp — sin q + Q sin q cos q, 

where q is the angle between pendulum and vertical, and 7 > is a friction coeffi- 
cient. This equation also describes the motion of a particle in a symmetric potential, 
shaped as a single well when < 1 and as a double well when Q > 1. The origin O 
is always an equilibrium, while for Q > 1, two new stable equilibria appear at 

Q± = (±g*(fi), 0), q*(to) = Arccos fi~ 2 . (47) 

The eigenvalues of the linearization of (pEBJ ) around O and Q± are given, respectively, 
by 

a° ± = -7 ± v / 7 2 + ^ 2 -l, < = -7 ± Vl 2 - fi2 + V- 2 . (48) 



16 



Figure 10. Solutions q(r) for slightly different values of the adiabatic parameter 
e. The time scale has been contracted in such a way as to show 20, resp. 40 periods 
of fi. One observes (a) solutions with the same period then f2(r), (b) solutions 
with twice the period of ^(r), going alternatively to one side and the other one, 
and (c), if e is carefully adjusted, solutions which have no apparent period. 



There are four qualitatively different phase portraits, delimited by the values Q = 1 
and Q = where 

(M 7 ) 2 = l-7 2 , fi + (7) 2 = H7 2 + V / 7+4], (49) 

namely (see Fig. 

• when < Q < Q^(j), O is a stable focus; 

• when ^-(7) < O < 1, O is a stable node; 

• when 1 < Q < ft + (j), O is a saddle and Q± are stable nodes; 

• when Q > ft + (j), O is a saddle and Q± are stable focuses. 

If O = fl(et) is made slowly and periodically time-dependent, we obtain the 
adiabatic system 



eq = p 

ep = —2jp — sin q + fi(r) 2 sin q cos q 



(50) 



This system displays two interesting phenomena. The first one is a bifurcation 
delay similar to the one already observed in previous examples: when Q is increased 
beyond 1, the pendulum remains for some time in unstable equilibrium close to the 
origin, before joining one of the stable equilibria Q + or Q_. When Q is decreased 
again below 1, the pendulum follows this equilibrium until it joins the origin again, 
leading to hysteresis. The second interesting phenomenon is related to the sequence 
of visited equilibria, which depends on the value of the adiabatic parameter (Fig. 



T0|) . For some values, the pendulum always chooses the same equilibrium, just as the 
Lorenz system always chooses the same direction of rotation for the convection rolls. 
For other values of e, however, one observes a sequence with twice the driving period, 
in which the pendulum visits alternatively the equilibria Q + and Q_. Between these 
periodic behaviours, it is even possible to observe apparently random sequences, 
which we called chaotic hysteresis |[BK|| . 

In order to explain this behaviour, we now compute an asymptotic expression 
for the Poincare map in the (g,p)-plane, during one period of Q(t). If fi(r) remains 
within the interval fi + ], the system can be reduced to ID as in Section H, and 
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Figure 11. Function f2(r) considered in the analysis. The instants when f2 crosses 
the values Q± and 1 delimit the different phases of the motion. 



there is no possibility for chaotic motion. We thus consider the case where f2(r) has 
a larger amplitude (Fig. |TTj). It is useful to introduce the notations 

«°(r 2 , n) = Re f 2 a° + (r) dr, 0°(r 2 , t x ) = Im P a° + (r) dr. (51) 

Similar functions a* and 0* are defined for the linearizations around Q±. 

For f < t < 1, orbits are attracted by the stable origin. They remain close to it 
until a bifurcation delay time f + 1, defined by the relation o°{t + 1, f) = 0. During 
this part of motion, the system can be essentially described by its linearization 
around the origin. Except in a neighborhood of t±, where the eigenvalues a°± cross, 
we can carry out a dynamic diagonalization as in Section |[ The actual crossings 
are described by a local analysis, using Airy's equation. Combining these steps, we 
obtain that 

c ^ *,* is( cos(^) e- 5 2°/ £ sin(^ + 6°)\ w /cn . 

where <p° = 0°(t^, t°_) + 0(e) is the dynamic phase of oscillations around the origin, 
and the columns of S(t) are close to the eigenvectors associated with the origin. 
The positive factors 8° describe the asymmetric contraction due to the difference 
between a° + and a°_, and the 8° are geometric phase shifts. It can be shown that the 
effect of nonlinear terms can be absorbed in these small geometric corrections. 

The part of motion between f and f is essentially nonlinear. Near the origin, we 
may use invariant manifolds as in Section ^| to transform (|50| ) into 

ei = [a° + (T,e)+{3 + (£,r},T,e)}£, 
er] = [a°_(r,e) + /3-{£,r],T,e)]ri, 

where a±(r,e) = cl±{t) + 0{e) and j3± are of order |£| + \r}\. Starting at f with a 
small initial condition (£ > 0, r) ), the second equation in ([53]) shows that rj becomes 
exponentially small. The first one is used to prove that £ reaches a distance d from 
the origin (d not too large) at a time r(£ ) + 0(e), where 

a°{f{ti ),t) = -eHt /d), (54) 
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Figure 12. Schematic shape of the function Ti(£) of equation (|56| ) (a) for 
cos(0°/e) > and (b) for cos(0°/e) < 0. In the first case, there are two sym- 
metric stable fixed points. In the second case, there is a stable orbit of period 
2. 



provided £o > £ c = e~ a °^ T +^l £ (for smaller £o, the orbit does not reach Q + before 
the time r*). For r > f(£ ), the trajectory is attracted by Q + , around which we 
carry out a similar analysis than around the origin, with the result 

£(f)=<7e 1/4 + <^ /e cos(£) 
77(f) =e (a *-' 5 * )/e sin(f + 0*), 

where a* = a*(f,f(£ )) + C(v^) and 0* = 0*(r|,f(£ o )) + 0(y/e), while C, 5* and 
9* are constant at lowest order in e. The position at f thus depends essentially on £0 
through the delayed bifurcation time f (£0), in such a way that (|55|) is the parametric 
equation of a squeezed spiral (which is essentially the image of the unstable manifold 
of the origin under the flow from f to f). 

Combining this result with (|52|), we finally obtain a Poincare map of the form 



6 = m ;vo,e) =cos(f ) [C^ 1 / 4 + e Q */ £ cos(f )" 

+ e^- 5 ^ sin(f + 9°) sin(f + 0*) , ( 56 ) 

77 1 = T 2 (£ ;77o,£) ^(e^/e), 

where ^ = £(f + j) measures the distance to the stable manifold of the origin (it 
is close to q), and r]j = r](f + j) measures the distance to the unstable manifold. 
This expression is valid for £ > 0, but is easily extended to negative £0, since the 
Poincare map is odd. The dynamics is thus essentially determined by the ID map 
£0 | — ► Ti(£ Q ;0,e), which is oscillating around ±Ce 1 ^ 4 cos(^) (Fig. |12|). One easily 
shows the existence of a positive constant fx such that, if e 1 ^ cos((f) /e) > e -M//e , this 
map admits stable fixed points at ±£* « ie 1 / 4 cos(0°/e), which correspond to cycles 
of period 1. When e 1 / 4 cos{4>° / e) < — e _M//e , there is an orbit of period 2, for which 
the pendulum alternatively visits the left and right equilibrium (Fig. |1^). These 
properties can be shown to hold for the 2D map, which is confirmed by numerical 
simulations (Fig. |T3| ). 

Chaotic motion is possible in the intermediate regions, where |e 1//4 cos(0°/£:)| < 
e~w £ . p or the 2D map, it is difficult to prove existence of such a motion, but one 
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Figure 13. (a) Numerically computed bifurcation diagrams of the Poincare map. 
For each value of e, we have plotted the asymptotic value of q(f + n), for one 
initial condition. On the domain < e < 0.025, the diagram clearly shows the 
alternates of regions with a period-1 and a period-2 cycle, separated by small 
chaotic zones, (b) Plots of the functions e 1 / 4 cos(0°/e) and ±e~ M / £ . Light gray 
zones are those where the theory predicts existence of a period-1 cycle, medium 
gray zones those with a period-2 cycle. Dark gray zone are those where chaotic 
hysteresis is possible, and, indeed, observed. We point out that in figure (b), the 
dynamic phase (f>°(0) has been computed analytically. Only the next-to- leading- 
order correction to 4>°(e) (which results in a phase shift) has been chosen in order 
to fit the numerical results. 




h 
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Figure 14. The interval map (a) admits the Markov subgraph (b), which allows 
for periodic orbits of all periods, except possibly 3. 
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can do more for the simplified ID map, using symbolic dynamics. In fact, when 
cos(4>°/e) = and under certain conditions on Q(t), one finds that Ti(£) behaves 
as in Fig. [I4|a: it vanishes at two points £1 and £2, and, being odd, also at — £1 and 
—£2- These points define four intervals I-2, I-i, h and I2, and the maximum of T\ 
on 1\ is larger than £2, while its minimum on I2 is smaller than — £1. The Markov 
graph of T\ is defined as the graph with sites Ij, admitting an oriented edge Ij — > Ik 
whenever T±(Ij) D Ik- It is known ||MS|| that for every path in the Markov graph, 
there exists an orbit visiting the corresponding sequence of intervals. In particular, 
using Sarkovskii's theorem [ |MS|| , it is possible to prove the existence of periodic 



orbits of every period except possibly 3. 

The mathematical pendulum analysed here has been realized experimentally, and 
all the phenomena predicted by the equations have been observed. They depend, in 
fact, only on a few qualitative features of the system. The origin should be a focus 
for some values of the parameter, in order to allow the orbits to wind around it. It 
should be hyperbolic for other values of the parameter, for which two new stable 
equilibria should exist. Chaotic motion requires, in addition, these asymmetric 
equilibria to be sometimes focuses, in order to create the oscillations in the Poincare 
map. Under these conditions, it should be possible to observe chaotic hysteresis for 
other nonlinear oscillators (see Section ISTi 



6 Examples of eigenvalue crossings 

We mentioned in Section [3] that adiabatic linear systems of the form ex = A(r)x 
could be diagonalized (and thus solved) by the change of variables x = S(t; e)y, 
where S is a matrix satisfying the equation 

eS = AS - SD, (57) 

and D is a suitable diagonal matrix. This procedure is only useful, however, if we 
manage to control the transformation matrix S, which should be bounded (e.g. close 
to the matrix So(r) which diagonalizes A statically). Such a control turns out to be 
possible at least in two cases: when the eigenvalues of A(r) have different real parts, 
or, (in a more restricted sense) when they have the same real part but different 
imaginary parts.0 

This leaves open the question of the effect of different types of eigenvalue cross- 
ings. The most generic case, when A is not diagonalizable at the crossing time, has 
been mentioned in Section ^. It can be studied using the properties of Airy func- 
tions (see also |Wa|| ) . In the present section, we illustrate the effect of two other 
types of crossing. The first one occurs when A{r) is symmetric, and can thus be 
diagonalized even when it has identical eigenvalues. The second one arises when the 
eigenvalues' real parts cross, but their imaginary parts are different. We call this 
situation eigenvalue cruising; it is closely related to properties of dynamic Hopf 



bifurcations discussed in [Nel, Ne2 



4 More generally, the system can be bloc-diagonalized when the eigenvalues can be split into 
two groups with non-crossing real parts. 
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6.1 Symmetric crossing 

Let us consider the overdamped motion of a particle in the 2D potential 

$(x,et) = -~(x\A(et)x) + i||x|| 4 , (58) 
where A is a symmetric matrix. The equation of motion can be written 

ex = A{t)x — ||a;|| 2 a;. (59) 
We assume that the matrix A(r) is given by 

A(r)-a(r)( COs26{T) sin2 ^A (Q ) 
A(r) - a[r) ^ _ cog 26{r) J , (oUj 

so that it admits eigenvalues ±a(r) and eigenvectors v\ = (cos 9, sin 9) and v 2 = 
(— sin#, cos#). Thus, the potential $ has minima at ±^/aVl if a is positive, and at 
±y/—a v 2 if a is negative. To diagonalize the linearized equation ex = A(r)x, we 
may try to solve equation fl5?D with matrices S and D of the form 

„_/cos0i(r) ~sm9 2 (r)\ _ (d x {r) \ . , 

\sm9 1 (r) cos0 2 (r);' U \ d 2 (r) J ' 1 ' 

Substitution of this Ansatz in fl5"T|) yields the relations 

e9 1 ^-a(r)sm2(9 1 -9(r)), d^r) = a(r) cos 2(9 1 - 0(t)), (62) 

e9 2 = a(r) sin 2(9 2 - 0(t)), d 2 (r) = -a(r) cos 2(fl 2 - 0(r)). (63) 

If a(r) does not vanish (i.e., when there is no eigenvalue crossing), these equations 
admit equilibrium branches at 9\ = 9 2 = 9(t), of opposite stability. By the results of 
Section^], we know that they admit particular adiabatic solutions 9\{r) = 9(t)+0(s) 
and 9 2 {j) = 9{r) + 0(e). The evolution operator of the linearized system can thus 
be written 

U(r,r ) = S(r)( e ^ )/£ ) S(r )-\ 

r V r J (64) 

5i, 2 (r,r )= / d 1>2 (s)ds = ± a(s) ds + 0(e 2 ). 

J TO J TO 

The columns of S(t) can be considered as dynamic eigenvectors which are close 
to the static eigenvectors v\ )2 . They define invariant subspaces (depending on r), in 
which the motion is expanding, resp. contracting. 

When a(r) is allowed to vanish, new phenomena occur because the equations 
(p3) undergo bifurcation. It is instructive to consider the case a(r) = — cost, 
for three different functions 9(t) : (1) 9(t) = 0, (2) 9{r) = — cost and (3) 9{r) = r. 

If 9(t) = 0, (p2]) admits the solution 9\ = 9 2 = 0, and d± t2 = ±a(r). The 
evolution operator can thus be written 

/ p -sin(r)/e n \ 

^(r,0)=( Q eS in W /J- (65) 
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Figure 15. Solutions (thin lines) of ( p^) when a(r) = — cost: (a) in the case 
9(t) = — cost, (b) in the case 9(t) = r. In both cases, one can construct a par- 
ticular solution remaining close to the static equilibrium 6{t) (thick lines, where 
the solid lines indicate stable branches and the dashed lines unstable ones), admit- 
ting a discontinuity at r = y. Solutions of ( |63| ) behave in a similar way, with a 
discontinuity at r = S. 



The subspaces X\ — and x<i = are invariant. If £2(0) = 0, xiij) remains 
exponentially small until r = 7r, which is the standard bifurcation delay. 

If 6{t) = —cosr, it is not possible to construct solutions of ( |62|) remaining 
indefinitely close to 9{r). The best one can do is to construct periodic solutions 
6\{t) and 6 2 (t) admitting a discontinuity of order y/e, respectively at times 4f and 



I (Fig. p~5|a) . As a result, for ~ < r < 4^ we have 



/ 5i(r,7r/2)/ £ q \ / 5 1 ( 7 r/2,0)/ £ n \ 

^,0)=5(r)^ Q e , 2(W2)/e J T ^ Q e52(7rA o )/£ J ^O)- 1 , (66) 
where 

r = s ( | + )-'5(i-) = J 1 + ° <^> - f + °M) , (67 ) 

with $2 = 6 I 2(|±). The off-diagonal term of this matrix induces a transition between 
the directions which were invariant before r = |. In particular, when r = 7r, we 
have to leading order in e 

^,o)«(J -<*-«.-)<"•). (68) 

This transformation rotates the vertical axis by almost tt/2. The matrix U(2tt,tt) 
is found to rotate the horizontal axis by almost — tt/2. This means that there 
exists no invariant subspace in which the particle performs an independent motion. 
The eigenvalue crossing thus results in an interaction between both modes, with the 
particle always following the most unstable direction. The sign of the discontinuities 
of 6*i 5 2 is important: in this case it induces a back-and-forth oscillation of the particle 
between two wells (Fig. |T6|a). 



23 



Vl 





V2 







Figure 16. Solutions of (|58| ) when a(r) = — cost, plotted with respect to the 
rotating reference frame (yi,v 2 ): (a) in the case 9(t) = — cost, the particle oscil- 
lates back and forth between two potential wells, while (b) in the case 6(t) = t it 
visits all four wells in a row. 



In the case 9(r) = r, the situation is 
opposite sign (Fig. |15|b). As a result, the 
same direction, and the particle visits all 



similar, but with a discontinuity of 6\ of 
coordinate axes are always rotated in the 
four wells in a row (Fig. p~6|b) . 



6.2 Coupled oscillators and eigenvalue cruising 

We call eigenvalue cruising the situation arising when some eigenvalues of the 
matrix A drift past one another at some imaginary distance. This cruising also 
leads to an interaction between the modes; however, unlike in the case of diagonal 
crossing, this interaction is delayed. 

Eigenvalue cruisings appear in particular in coupled oscillators. Consider for 
instance the system 

q\ + 2 7 ift + (1 - A + q\ + Aq 2 2 ) qi - /iq 2 = 0, 
q 2 + 2 l2 q 2 + 4(4 - A + q\ + Aq 2 2 )q 2 + m = 0, 



which was introduced by Kobayashi [ |Ko| to describe the vibrations of a buckled 



plate with supersonic flow on one side of the plate. The variables q\ and q 2 are 
amplitudes of the two dominant Fourier modes of the deflection, A is the in-plane 
compressive stress, \x the dynamic fluid pressure of the supersonic flow, and 71,2 are 
friction coefficients (which were taken equal in [|Ko| ) . 

Introducing pi = qi and p 2 = q 2 , ( |69| ) can be written as a 4D first order system 
for the variables {qi,Pi,q 2 ,p 2 ), which admits the origin as an equilibrium. The 
linearization around the origin is a 4 x 4 matrix with eigenvalues 



a 1)± = -7i± V7i 2 + A-l + CV), 

(70) 



a 2>± = - 72 ± ^ 72 2 + 4(A - 4) + 0(fi 2 ). 

An eigenvalue cruising arises for instance in the following situation: assume 71 = 2, 
72 = 1 and fi = 0, so that a\^± = — 2 ± VA + 3 and a 2 ,± = — 1 ± i \/15 — 4A. As 
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Figure 17. Eigenvalues of the linearized Kobayashi equations for 71 
and [i = 0. There is an eigenvalue cruising for A = — 2. 



2, 72 = 1 



A increases from — 3 to ^f, the complex eigenvalues a 2j ± correspond to oscillations, 
while the real eigenvalues a%,± describe an overdamped motion. There is a cruising 
at A = —2 and the origin becomes unstable at A = 1 (Fig. [IT]). The same qualitative 
features hold for small positive coupling \i. 

We are interested in the following question. Assume that A is increased mono- 
tonically and adiabatically, starting with a generic initial condition at a time tq 
where A is smaller than 1. For what value of A does the trajectory depart from the 
origin? The answer turns out to be related in a rather subtle way to bifurcation 
delay and eigenvalue cruising. It is easier to explain this phenomenon on the simple 
model equation 

As in Kobayashi's equations, a\ represents the overdamped mode, and a 2 the oscil- 
lating one (in complex notation). The cruising occurs at r = 0. To diagonalize this 
equation, we try to solve the equation eS = AS — SD with the Ansatz 

*M = U) 1 T) ) • W = *w) ' (72 » 

Substitution in the equation for S yields the relations 

esi = —ji — (r — i)s\ — fxsl, d\ = at + fisi, (73) 

es 2 = [i + (r - i)s 2 + ns\, d 2 = a 2 - /xs 2 . (74) 

The first equation has a static equilibrium at s*(r) = — /i/(r — z) + (9(/x 2 ), which 
is unstable for r < and stable for r > 0. One can show that the solution of (|73| ) 
with initial condition si(0) = s*(0) tracks the branch s*(r) at a distance at most 
(9(e), for both negative and positive times. 

The second equation has a more subtle behaviour. It admits an equilibrium 
branch at s 2 (t) = /x/(r — i) + 0(fi 2 ), which is stable for r < and unstable for 
r > 0; in fact, it undergoes Hopf bifurcation. Such bifurcations have been studied 
by Neishtadt [Mel , [Ne2| , [Ne3|| . The interesting fact is that there exist solutions 



tracking the equilibrium branch beyond the bifurcation point, but only until a time 
called maximal delay or buffer point. This point is obtained in the following 
way: let a(r) = r — i+0(fi 2 ) be the linearization of (fF4[) around the equilibrium 
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Figure 18. (a) The level lines of the function Rea(r) of equation (|7^) are hy- 
perbolas centered at r = i. The largest positive time which can be connected to 
the negative real axis by such a line is the buffer time = 1. (b) The origin 
becomes unstable when the largest exponent of a matrix element of the evolution 
operator (77) becomes positive. This may happen earlier when [i ^ 0, because the 
oscillators interact at the buffer time. 



Sofr). Define the function 



Rea(r) = Re/" a(s) ds = \ [(Rer) 2 - (Imr - l) 2 + l] + £>(/i 2 ). (75) 
Jo 

The buffer time Tb is the largest real time which can be connected to the negative 
real axis by a path of constant Re a (and with some additional properties given in 
Ne3|| ). In the present situation, — 1 (Fig. |18|a). As a consequence, for r ^ 1, we 



can construct a solution of ([7|) which is close to ^(t), an d the evolution operator 
of ( |7TD is given by 



U{t, 


ro) 


= S(t) 








Si(r, 


to) 


















ro) 


= f d 






J T 



e 5i(T,ro)/e Q 

Q e S2(r,r )/e 



d 1 (s) ds = |(r 2 - r 2 ) - (r - r ) + 0(fi 2 ), (76) 

d 2 (s)ds = (i-l)(r-r ) + O( l x 2 ). 

For r > 1 > tq, however, the solution S2(t) necessarily admits a discontinuity of 
order e/i at the buffer time. A similar calculation as in the previous subsection yields 
the evolution operator 

/ p h(r,T )/e r)( FII \JS 1 (T,l)+S 2 (l,T )]/e\ 

U(t,t ) = S(r) 1° Q [ ^ e s 2{ T,T 0)/ e JSiro)- 1 (77) 

If the system starts away from the origin at ro < 0, it will follow the origin exponen- 
tially closely until a delay time f, which is the first time at which one of the matrix 
elements becomes of order 1 again. If fi = 0, we simply have f = 2 — r . When 
7^ 0, however, we cannot overlook the interaction between the overdamped and 
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the oscillating mode, which takes effectively place at the buffer time = 1, and 
may cause the system to become unstable at an earlier time. For instance, when 
r = —2, the usual delay time for fi = would be f = 4, while the effective delay 
time for p > is 1 + V§ + 0(/i 2 ) (Fig. 0b). 

A similar phenomenon is observed for the Kobayashi equations (0), only with 
different values of the cruising, buffer and delay times. The effective value of the 
bifurcation delay can be shown to depend at leading order in e only on the lineariza- 
tion around the origin. It is important, together with nonlinear terms, for the global 
structure of motion, since it influences the choice of the asymmetric equilibrium the 
system follows after leaving the origin. In fact, for large amplitude oscillations of 
the form A(r) = 8sin(r), we observed numerically that the Kobayashi equations 
display chaotic hysteresis just as the rotating pendulum in section |5|. This is not 
really surprising, since even when fi = 0, each oscillator is similar to the rotating 
pendulum of the amplitude of A is large enough. A positive fi, however, will modify 
the bifurcation delay and the dynamic phases and amplitudes which determine the 
structure of the Poincare map. 
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